Chasing the non-hnear evolution of matter power spectrum with numerical 
resummation method: solution of closure equations 



Takashi Hiramatsu^ and Atsushi Taruya^^'^ 

^Institute for Cosmic Ray Research, The University of Tokyo, Kashiwa, Chiba 277-8582, Japan 

^Research Center for the Early Universe, School of Science, 
The University of Tokyo, Bunkyo-ku, Tokyo 113-0033, Japan and 

^\ ' ^Institute for the Physics and Mathematics of the Universe, 

' The University of Tokyo, Kashiwa, Chiba 277-8568, Japan 

o ■ 

CNj ■ We present a new numerical scheme to treat the non-hnear evolution of cosmological power 

spectra. Governing equations for matter power spectra have been previously derived by a non- 
', perturbative technique with closure approximation. Solutions of the resultant closure equations just 

correspond to the resummation of an infinite class of perturbation corrections, and they consistently 
reproduce the one-loop results of standard perturbation theory. We develop a numerical algorithm 
to solve closure evolutions in both perturbative and non-perturbative regimes. The present nu- 
^is^ . merical scheme is particularly suited for examining non-linear matter power spectrum in general 

cosmological models, including modified theory of gravity. As a demonstration, we study weakly 
non-linear evolution of power spectrum in a class of modified gravity models, as well as various dark 
energy models. 
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C/3 , In the last decade, the late-time cosmic acceleration has been one of the most important discoveries in physics 
, ^/ and cosmology (e.g., Refs. P, Q)- Although the origin of late-time acceleration is thought to be a mysterious energy 
component called dark energy, a possibility of long-distance modification of general relativity is still viable (e.g., 
(M ■ Refs. Hi]), and our understanding of the nature of late-time acceleration is still lacking. Currently, the dark energy 
^ ' equation-of-state parameter Wdc, which is phenomenologically introduced to characterize the cosmic acceleration and 
is defined as the ratio of the pressure to the energy density of dark energy, is consistent with a cosmological constant 
^T" ' (''"do = —1) at a level of 10% precision, and with no evidence for time dependence (e.g., Refs. 6]). Toward a deeper 
, understanding of the nature of late-time acceleration, a precise measurement of both the cosmic expansion history 
and the growth of cosmic structure is a key to to distinguish between different models of dark energy, as well as to 
discriminate the dark energy from the modification of gravity. 

Among several observational techniques, baryon acoustic oscillations imprinted on matter power spectrum and 
I cosmic shear, measured from galaxy samples, are the most promising techniques sensitive to the expansion history 
and growth of structure. A crucial remark is that they strongly rely on the accurate predictions of non-linear 
matter power spectrum. Hence, in addition to the precise measurement, a high-precision theoretical template for 
the non-linear power spectrum must be developed in order to achieve order-of-magnitude improvement of the current 
^ ' constraints. 

5^ [ Recently, several analytical approaches to predict the non-linear power spectrum have been developed, comple- 
mentary to the N-body simulations [Z, M, H flO - fill . [H, [Hi] . In contrast to the standard analytical calculation with 
perturbation theory (for a review, see Ref. these have been formulated in a non-perturbative way with tech- 

niques resumming a class of infinite series of higher-order corrections in perturbative calculation. Thanks to its 
non-pcrturbativc formulation, the applicable range of the prediction has been greatly improved, and the non-linear 
evolution of baryon acoustic oscillations was found to be accurately described with a percent-level precision. 

Note, however, that these analytical calculations involve several approximations or simplifications in order to make 
the analysis tractable. This severely limits the applicable range and/or the versatility of predictions. For example, 
in Refs. [12,1^, a perturbative treatment called Born approximation has been partly adopted in order to evaluate 
the non-perturbative expressions for power spectrum. Furthermore, most of the analysis presented so far rely on the 
Einstein-de Sitter approximation, in which all the calculations done in the Einstcin-de Sitter universe are extended to 
apply to the other cosmological model by simply replacing the linear growth factor in Einstein-de Sitter universe with 
that in the other cosmology (see Sec. IV B II in detail). This is very crucial in studying the non-linear matter power 
spectrum in general cosmological models, especially in modified gravity models. 
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In the present paper, in order to bring out the advantage of non-perturbative formulation as much as possible, we 
present a numerical resummation scheme to calculate the non-linear matter power spectrum. Our treatment relies on 
the formalism developed by Ref. [13], in which the non-linear statistical method used in the subject of turbulence (e.g., 
Ref. [l5l|) was applied to the derivation of governing equations for power spectrum. The resultant equations called 
closure equations are the non-linear intcgro-differential equations coupled with non-linear propagator. The solution 
of closure equations effectively contains the information of the higher-order corrections, similar to the renormalized 
perturbation theory by Crocce & Scoccimarro 0, M, 113 • 1^ been shown that the analytical predictions based on 
the leading-order Born approximation agree with N-body simulations very well in a mildly non-linear regime, and 
a percent-level precision was achieved at some ranges 11611 . The agreement of the prediction is further improved if 
taking account of the next-to-leading order correction [l7[. Hence, with the numerical implementation of the closure 
equations, all orders of Born approximation are included, and the prediction will be much better than the analytical 
treatment. Further, the numerical treatment is particularly suited for studying the non-linear power spectrum in 
various cosmologies where the analytical calculations with Einstein-de Sitter approximation is no longer possible. 

The paper is organized as follows. In Sees. Ull and [TTTl we briefly review the basic treatment of our approach and 
formalism. We then discuss how to solve closure equations in Sec. IIVI As shown in Ref. , the closure equations 
automatically reproduce the leading-order results of standard perturbation theory if replacing the quantities in non- 
linear terms with linear-order ones. This treatment has been used for computing quasi non-linear spectrum in modified 
gravity models in Ref. [l^. In Sec. |Vl we present numerical solutions of closure equations in both full non- linear 
and perturbative treatment and demonstrate how the present scheme can treat analytically intractable cases. Finally, 
Sec. I VII is devoted to the discussion and conclusion. 



II. PRELIMINARIES 



Throughout the paper, we consider the evolution of mass distribution in the flat universe, neglecting the tiny 
contributions from the massive neutrinos. We treat the cold dark matter (CDM) plus baryon system as a pressureless 
perfect fluid. Then, assuming the irrotationality of fluid flow, the governing equations for matter distribution become 
the continuity equation and the Euler equation coupled to the Newton potential (p (e.g., Ref. : 

dS(T,x) , ^ 1 . 



dr \ dT J ' ' ' a?m ' ' 

where 6 is the mass density field, and 9 is the velocity divergence defined as = V • v/{aH). Here, we introduce the 
time variable given by r = log(a/ao), with qq being the scale factor at the present time. With this time variable, the 
fiat Friedmann equation becomes 



= Hi |r!,ne-3- + ride exp[-3^ (It' {I + wa,{t')}\ | . (3) 



The quantity Hq is the Hubble parameter at the present time, and and Side are the density parameters of the 
matter and dark energy, respectively. 

To treat the non-linear evolution of matter power spectrum, we will work with the Fourier transform of the fiuid 
equations, ([1]) and They are given by 

^'^'^'"^ ■ ^^(fc,r) = - /^^!^i^Mfc-fci-M(l + ^^(fc^^ (4) 



dr ' ^ ' ^ J (27r)3 ' \ ' |fci|2 



d9{k,T) ( d\nH\ ( k 
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J - fci - k2) 9ik„r)9ik2, r). (5) 



As for the Poisson equation, we have 



^{k, t) = 4^ Geff (fc, t)p„, 5{k, t). (6) 
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Here, GcS is the effective Newton constant, which genericaUy depends on the scale and time in modified theory of 
gravity. In principle, the Newton potential can be a non-linear function of the density field. In fact, successful modified 
gravity models that explain late-time acceleration such as the Dvali-Gabadadze-Porrati (DGP) model 3 and f{R) 
gravity models (for a review, see Ref. have non-linear interaction terms, which are essential to recover the general 
relativity on small scales |l9j . In the present paper, we restrict our analysis to the cases with linear Poisson equation. 
The extension to the non-linear case is straightforward and is discussed in a separate paper [l8l |. 

The evolution equations ([3]), (O and ^ can be further reduced to a compact form if we introduce the following 
quantity: 

$„(fc,T)= (5(fc,r), -0(fe,T)); (a =1,2). (7) 
Then, we write down the evolution equations in a single form as 

where ^acd is the vertex function defined as 



Aa6$b(T, k) = II ,^l,g fo(fci + fc2 - fe)7acd(fcl, fe2)*c(fel)$d(fc2), (8) 



7112(^2, fci) = 7i2i(A;i, A;2) = ^ (l ' 



7222(^1, ^2) = 



2 V IfeiP 
1 f\ki + k2\'^ki ■k 2 

|fclP|fc2P 



(9) 



The operator Kab is defined by 



with the matrix fJ^b being 



Aab = Sab^ + ^abik,T), (10) 
OT 



^ab{k, r) = i^^Q^^Prn ^ , d\n H ) . (11) 




III. CLOSURE EQUATIONS 



In the present paper, we are especially concerned with the non-linear evolution of power spectrum defined by 

($a(fe,r)$6(fc',T)) = i27Tf5D{k + k')Pab{\k\:T), (12) 

where the bracket (•) stands for ensemble average. In the above definition, we have the three different power spectra, 
Pii, P12 = P21, and P22, which respectively correspond to Pss, —Pse and Pgg. 

For the analytical calculation of the power spectrum, the standard treatment of perturbation theory is to expand 
the quantity as $q = <I>i^^ + + • • ■ , and to iteratively obtain the solutions from Eq. ([8]). Substituting 
the perturbative solutions into the definition (|12p . we obtain the non-linear corrections to the power spectrum. This 
treatment is straightforward, but is plagued by a poor convergence of the perturbative expansion. Because of this, 
the applicable range of the standard perturbation theory (SPT) is restricted to a narrow range on large scales. 

Recently, the improved treatment of the perturbation theory has been proposed by several authors employing the 
so-called renormalized/resummation techniques 0, H, H, [13, [U, Il2j IH, US, EU, H^l- In these treatments, the naive 
expansion of the SPT is re-organized by introducing the non-perturbative statistical quantities, and the information 
of the higher-order corrections in SPT is effectively incorporated into each order of expansions. As a result, even 
truncating the expansion at some orders still contains the non-perturbative effects of non-linear clustering, leading to 
the improvement of the convergence properties. 

Here, among several non-perturbative techniques, we consider the closure theory proposed by Ref. in which 
we have applied the non-linear statistical method commonly used in the subject of turbulence (e.g., Ref. [ill) to the 
cosmological perturbation theory. In this treatment, the renormalized expansion has been first constructed according 
to the renormalized perturbation theory by Ref. [3- Then, we truncate the expansions at the one- loop order. Under 
the tree-level approximation of the vertex function, this leads to a closed system of the power spectrum and non- 
linear propagator. Though some non-perturbative properties are missed in this treatment, an advantage of this 
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formulation is that we can compute the power spectrum numerically by solving the evolution equations, keeping 
full non-perturbative information of the non- linear clustering at the one- loop order. This forward modeling may be 
suitable for a fast computation of the power spectrum, unlike the backward treatment of the perturbative expansions, 
which requires the time-consuming multi-dimensional integrations. 

Let us define the non-linear propagator, Gab(A;|T, t'), and the cross power spectra between different times, 
Rab{k;T, t'): 

S<^a{k, t) ^ ^ Gab{k\r. r')SDik - k'), (13) 



($,(fe,T)$f,(fe',T')> = (2^)3<5i5(fc + fe')^a6(|fc|;r,r'); (r > r'). (14) 

Then, the governing equations for power spectrum, equivalent to the rcnormalized expansions truncated at the one- 
loop level, become [IS] 



^abcdPcd{k]T) ^ / dr" MUk;T,T")Rbs{k-T,T") 

J To 

+ I dr" Na,{k-T,T")GM{k\T,T") 



+ (a ^ 6), (15) 
AabRbcik;r,T')=. f dr" Mas{k;T,T")Rsc{k;T",T') 



+ / dT"Nai{k;T,T")G,iik\T',T"), (16) 

J To 

AabGbc{k\T,T')= f dr" Mas{k;T,T")Gsc{k\T",T'). (17) 



Here, i?^(fc;T",T') = Rsc{k;T" ,t') for r" > r', Rcs{k;T' ,t") for r" < r'. The operator Sa&cd is defined by 

d 

^abcdir) = SacSbd^ + 5aAd{T) + Sbd^acir). (18) 

The matrices Mat and Nab are 

f d^k' I , 

Mas{k]T,T") = A J j^-^-fapq{k - k' ,k')jirsik' - k,k) 

X Ggi{k'\T,T")Rpr{\k-k'\;T,T"), (19) 

Naeik;T,T") = 2 J j^-^japqik - k' ,k')jirs{k ~ k' ,k') 

X Rgsik';T,T")Rpr{\k~k'\;T,T"). (20) 

Note that we have recast the original equations in Ref. [isj in more symmetrical way by changing the integration 
variable [c.f. Eqs. (49)-(53) of Ref . [H]]. By definition, the non-linear propagator and the cross power spectra should 
satisfy the boundary condition: 

Gabik\T,T) ^ dab, (21) 

lim Rabik;T,T') ^ Pab{k;T). (22) 

T — yT 

The closure equations p5p - (|17p are the integro-differential equations involving several non-linear terms, in which 
the information of the higher-order correction in SPT is encoded. Thus, replacing the statistical quantities Rab 
and Gab in non-linear terms with linear-order ones, the solutions of closure equations automatically reproduce the 
leading-order results of SPT, i.e., one-loop power spectra. Here, the linear-order quantities denoted by Rai, and G^j, 
satisfy 

A,bG^^(fc|T,T') = 0, (23) 

AabRt{k;T,T') ^0; (t>t'). (24) 

For the rest of this paper, we focus on the numerical treatment of the closure equations and demonstrate the evolution 
of matter power spectrum in both non-linear and quasi-linear regimes by changing the treatment of non-linear terms. 
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IV. NUMERICAL METHOD 



The closure equations (|15[) - (|17p are the non-linear coupled equations involving the time-consuming integrals over 
space and time. In order to numerically treat these messy integrals while keeping computational cost, we implement 
the method used by Ref. [9|, in which the propagator and power spectra are expanded by a set of basis functions of k, 
and integrated with respect to k in advance of the time evolution. We adopt the trapezoidal rule for the integration 
with respect to t and fc, and the central difference formula for the time evolution. To be precise, we first prepare a 
discretised set of k labeled as km for m = 1, • • • , Nk, where we denote fci and fcAr^. by /cmin and A:niax, respectively. We 
define a set of triangular-shaped functions as the basis functions: 

k kjn—i ^ u u 

Tm{k) = I krn+l - k <k<k (25) 

0; otherwise. 
Then we expand the non-linear propagator, the auto- and cross-power spectra as 

Pab{k';T) ^J2'P-b,m{T)Tm{k'), (26) 

m 

Rab{k';T',T) = ^7^ah,™(T',r)r^(fc'), (27) 

771 

Gabik'y, r) = ^ GabM^', T)Tm{k'). (28) 

m 

The above expressions together with basis function (|25|) imply that the power spectra and the propagator between 
the discrete points are evaluated by the linear interpolation according to the definition of the basis functions. Note 
that these functions do not satisfy the orthogonality in the sense that the integration of the product 7^(fc)7^(fc) over 
the continuous space of k does not vanish even if m 7^ n = m + 1. The set of Tmik) has the orthogonality only on the 
discrete space because Tmih) = S„ii is satisfied. 

Substituting Eqs. ([26|) - ((28)) into Eqs. and ((20)) . we obtain a separable form of the matrices Mab and Nab '■ 



Mas{k;T,T") = J2 Sg,Mr,rnnprAr,r'')'r}p^L,n^Jk), (29) 

all indices 

7V,,(fc;r,T")= '^^s,„r{Ty')n,rAry')r^Zrs,rn,nikl (30) 



all indices 

where T^*^^ and T*-^-* are given by 



(Pk' 

(27r)3 

xr™(fc')r„(fc-fc'), (31) 

TZ\rs,m,n{k)=2j (fc - fc' , fc')7... (^ - fc' , fc') 

xTm{k')%,{k~k'). (32) 

The details of the description on the integrations PT|) and can be found in Appendix [Xj Although Eqs. (IFD) 
and (|32p seem to have many components, most of them vanishes because the vertex function 'jabc has only three 
non- vanishing components given in Eq. The relevant components in the summation are listed in Table [D 

Based on the essential points of our numerical treatment described above, we now consider how to solve the closure 
equations in a cosmological setup. Basically, we perform the following steps (see also Fig.[T]): 

1. Set T = Tinit (or z = Zinit) for sufficiently small value of Tinit < 0, where the universe is well-described by the 
Einstein-de Sitter (EdS) model, and impose the initial conditions: 

-Pah (fc; Tinit) = ^ab(fc; Tinit, Tinit) ^ ( j j ) -^""1(^)7 Gafc(A:|ri„it, Ti„it) = 5ab, (33) 
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(a, s) 


{apq,£rs) for M 


(1,1) 


(112,121), (121,121) 


(1,2) 


(112,112), (112,222), (121,112), (121,222) 


(2,1) 


(222,121) 


(2,2) 


(222,112), (222,222) 


{a,e) 


[apq,lrs) for A'' 


(1,1) 


(112,112), (112,121), (121,112), (121,121) 


(1,2) 


(112,222), (121,222) 


(2,1) 


(222,112), (222,121) 


(2,2) 


(222,222) 



TABLE I: List of non- vanishing components of the function T^pg Ws m ri(^) ^qs. (|29p and pO|l . 



where Pinit(fc) is the hnear power spectrum given at the initial time ri„it. Note that, in order to ensure the 
vahdity of this prescription, the initial condition should be imposed early enough so that the influences of the 
non-linearity and the transient from the initial condition can be neglected. Appropriate value of the initial time 
has been chosen based on the convergence test in Appendix [Bl 

2. Perform the integration over k for all components of T^pq ^^g „ „ and T^^J^^^ ^ ^ using the trapezoidal rule [see 



Eqs. dni]) and 

3. Suppose that the coefRcients Gab.p{Tm,Ti) and Tiab,p{Tm,Ti) for < i < to have been already obtained (filled 
circles in Fig. [T]), we perform the summation in Eqs. (|29|) and (pO| . and compute the kernels Mab{k; Tm, n) and 

Nab{k:,Tjn,Ti). 

4. Calculate the non-linear terms in the right-hand side of Eqs. ([T5)) - PT|) . We use the trapezoidal rule for the 
integration over time r". 

5. Advance the time step from t„i to Tm+i, and obtain a new set of arrays Gab{k\Tm+i,Ti), Rab{k]Tm+i,Ti) and 
Pab{k;Tm+i) for < i < m, applying the central difference formula to Eqs. p^ - p?)) (open circles in Fig. [T]). 
For the edge of arrays (shaded circle in Fig. [T]), the boundary conditions pT|) and ([22|) are used to obtain 

Gab{k\Tm+l,T,n+l) and Rab{k;Tm+l,T„i+i). 

6. Repeat the steps 3-5 until the time t„_|_i reaches the final time. 

Note that the trapezoidal rule at step 4 and finite difference scheme at steps 5 are explicitly written as follows. For 
Eq. (fT7)) with (r, r') = (Tp,Tg), we have 

Gab,riTp+l,Tq) - gab^r{Tp-i,Tq) I \r ( \ 

— h ^lab\Tp)ybc,r\Tp,Tq) 

Zl\t 

V (34) 

^ ^ ^ra^as{kr J ^p^ sc^ri^ra^ ^q) ^ 
m—q 

for p > q. In cases with p — q, the differentiation in left-hand side of the above equation is replaced with the first-order 
difference. Here, we set Cm = 1/2 for m = otherwise c„,. — 1. Eqs. (|15p and (|16p can be also written similarly as 
above. 

The above procedure can be also used for the calculation of the one-loop spectra in SPT. As we mentioned in 
Sec. Illli the calculation in the SPT additionally needs the solutions Gj^^ and given in Eqs. and ([M)) in 
advance. They are used at the step 3 to compute the integration kernels, Af^^ and N^^i,, which are evaluated from 
Eqs. (|19p and (|^D|) just replacing the integrands with linear-order ones. In the same manner, at the step 4, we replace 
Gab and Rab in the non-linear interaction terms with linear-order quantities, Gj^^ and R^i,, respectively. 



V. DEMONSTRATIONS 



In what follows, we present the results of numerical integration of closure equations. We first demonstrate the full 
non-linear calculation and present the results in Sec. IV Al Then, we move to discuss the perturbative treatment and 
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Ti= Tm+1 
Ti — Tm 



known 
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"Id 


o 






'Tm 


Tm+1 



Pah : time evolution 
Gab : boundary condition 



time evolution 



FIG. 1: Schematic draw of the numerical procedure to solve closure equations. Each circle represents both Qab,-p and TZab,p 
{Vab.p for Ti — Tm) OYi the discrcte time grid, just omitting the wavenumber dependence. The horizontal and vertical axes 
correspond to the second and third arguments of those quantities. In the steps 3 and 4 mentioned in the main text, we compute 
the matrices Mab{k; Tm, Ti) and Nabik; Tm, Ti) for < « < m, and evaluate the non-linear terms, depicted as filled circles at the 
right edge of the 'known' region. In next step 5, we advance the time step and obtain Qab,p{Tm+i,Ti) and TZab,p{Tm+i,Ti) for 
<i <m. (blank circles), and 'Pab,p and Gab,p for (rm+i,Tf,m + 1) (shaded circle). 



examine the weakly non-linear evolution of the power spectrum in dark energy and modified gravity models in Sec. lVBl 
The initial power spectrum Pinit(fc) is calculated from the linear transfer function in the flat ACDM model. We adopt 
the cosmological parameters determined from WMAP five-year results j6ij: A^(fco = 0.002Mpc~^) = 2.457 x 10"^ 
Us — 0.960, flni = 0.279, h = 0.701, for the amplitude of curvature perturbation, scalar spectral index, density 
parameter of matter, and Hubble parameter, respectively. Unless otherwise stated, we assume the dark energy with 
equation-of-state parameter Wdc = — 1- 

The parameters of our numerical calculations include the initial redshift Zi^a, the cutoff wave number fcmax, the 
number of time steps Nr, and the number of Fourier mesh Nk- We set Nt = 172 and = 200 with constant interval 
in linear and logarithmic scales, respectively. For the initial redshift and cutoff wavenumber, based on the convergence 
test in Appendix |B1 we chose Zinit = 200 and fcmax ~ 5/iMpc^^. With this choice, the numerical errors in the SPT 
calculation are reduced to a sub-percent level. 



A. Full non-linear calculation 



In the full non- linear treatment, the solutions of auto- and cross-power spectra as well as the non-linear propagator 
are simultaneously obtained from the closure equations at each time step. Here, for illustrative purpose, we first show 
the non-linear propagators, which clearly manifest the non-perturbative property of non-linear clustering incorporated 
into our formalism. 

Fig.[2]plots the non-linear propagator G'ii(fc|z, Zinit) as function of wavenumber given at different redshifts, z = 0.5, 
2 and 5 (from left to right). Clearly, the numerical results depicted as solid lines exhibit the damping oscillation, and 
asymptotically approach zero at fc ^ oo. The characteristic scale of the damping is shifted to a lower k for decreasing 
the redshift. These behaviors are marked contrast with the linear theory prediction depicted as dotted line. Note 
that the results including the leading-order correction (one-loop SPT) slightly improves the low-fc behavior, but they 
eventually become negative and diverge at fc ^ oo. In this respect, the damping properties seen in the numerical 
results can be regarded as the non-perturbative effect, which results from the resummation of infinite series of higher- 
order corrections. Indeed, the damping behavior in the non- linear propagators has been already confirmed in the 
N-body simulations [1, [2^ , and is essential for the accurate prediction of power spectrum [l^ . 

In Fig. [21 the dashed lines indicate the analytic results obtained from Ref . [l^ . Basically, these are the approximate 
solutions of Eq. ([T7]) constructed by matching the asymptotic solutions in the low-A: and high-Zc limits. Although 
the analytic results at lower redshifts slightly deviate from the numerical solutions, the overall agreement between 
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0.01 0.1 1 

k [h/Mpc] 

FIG. 2: The non-linear propagators Gii(fc|2, zinit) at z — 0.5, 2 and 5 from left to right. The vertical axis is normalised by the 
linear growth rate, D{z). The solid lines represent the numerical solution of the closure equations (|15|l -^ (|17p . The dashed line 
are the approximate solution derived in Ref. [13j , and the dotted line indicates the linear theory prediction. 




k [h/Mpc] k [h/Mpc] k [h/Mpc] 



FIG. 3: Time evolution of matter power spectra, Pn (fc), evaluated at 2 = 3, 1 and 0.5 from left to right panels. The vertical axis 
is normalised by the linear growth rate, D{z). The solid lines represent the numerical solution of the equations (|15p - (|17[) . The 
dashed and dotted lines are the analytic results including the corrections up to the first- and second-order Born approximation, 
respectively [l^. The squares with errorbar indicate the N-body simulations taken from Ref. [l7| ]. 

these two curves is remarkable. This may be an independent check for the stabihty of our numerical scheme, and the 
accuracy of our code seems comparable to or even better than the analytic calculations. 
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Now, in Fig. [3l we show the redshift evolution of the matter power spectrum, Pii(fc; z), obtained from the closure 
equations. For comparison, we also plot the N-body results taken from Ref. [13] ■ Solid lines represents the numerical 
results of closure equations, dashed and dotted lines are the results of analytic calculations including up to the leading- 
order and next-to-leading order perturbative corrections, respectively. Here, the analytic results were obtained based 
on the integral solutions of the closure equations presented in Ref. [l^ We employ the Born approximation 
to evaluate the integral solutions perturbatively. Although the analytical treatment is found to accurately describe 
the non-linear evolution of baryon acoustic oscillations with a precision of sub-percent level [TgI . [iTj , because of the 
perturbative treatment, applicable range of the analytic treatment is limited to a narrow range. As clearly shown 
in Fig. [3l the resultant power spectra rapidly fall off at some higher wavenumbers. By contrast, the power spectra 
obtained from the numerical calculation first trace the analytical results on large scales, and they extend over small 
scales without a sharp drop of the amplitude. Remarkably, the numerical results quite resemble the N-body results 
at < 1 (0.6)/i Mpc^^ for z = 3 {z = 0.5), and the agreement between these two results reaches the accuracy of 
~ 4%(8%) level. This is a clear manifestation of the fact that full non-linear treatment of the closure equations is 
indeed a non-perturbative way of calculating the power spectrum beyond the weakly non-linear regime. Hopefully, 
it would be a fast computational tool complementary to the N-body simulations. To clarify the usefulness of this 
approach, a more quantitative comparison between N-body simulations and our numerical treatment is needed. We 
will discuss this issue in a future work. 



B. Perturbative calculation 



In this subsection, we turn to focus on the perturbative treatment of the closure equations, by which all the 
quantities in non-linear terms are replaced with the linear-order ones. As we mentioned, this treatment automatically 
reproduces the one-loop results of SPT. Owing to the numerical treatment, we can address weakly non-linear evolution 
even when the analytical calculations are no longer possible. In Sec. IV B 11 we discuss the one-loop power spectra in 
dark energy models, and address the validity of the analytical treatment based on the Einstein-de Sitter approximation. 
In Sec. IVB 2[ we examine a class of modified gravity models with linear Poisson equation, where the effective Newton 
constant manifestly depends on scale. We demonstrate how the modification of the gravitational-force law affects the 
power spectra in weakly non-linear regime. 



1. Dark energy models 



The one-loop SPT has recently attracted renewed interest for an accurate modeling of large-scale structure. In 
particular, a precise measurement of baryon acoustic oscillations made by ongoing and/or upcoming galaxy surveys 
to probe the nature of late-time cosmic acceleration provide a strong motivation to use the one-loop SPT for an 
accurate template of matter power spectrum (e.g., Refs. [H, [H, |2l, [2^ ) . In these experiments, the required accuracy 
for theoretical template reaches at a percent level. 

In the analytic treatment of one-loop power spectra, the Einstein-de Sitter (EdS) approximation has been frequently 
used in the literature (e.g., Ref. and references therein). Under the approximation, the higher-order solutions 
of perturbation are approximately described by the linear growth factor D{z), and the resultant power spectra are 
schematically expressed as 



Patik; z) = D\z)Pl,{k) + D\z)P'J°°^{k) 



(35) 



Note that for dark energy models in general relativity, the EdS approximation is mathematically equivalent to solving 
the closure equations just replacing the matrix flab in the operators Aab and 'Sabcd with 



/ 







-1 \ 

din/ 
dr / 



(36) 



with the function / defined by / = dlnD/dr. 

Here, we consider two specific examples of dark energy models characterized by the equation-of-state parameter 
Wdc as [13, HI 



Wdo(a) = Wo + Wail - a), 



(37) 
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FIG. 4: The density (left) and velocity divergence power spectra (right) normalized by the linear theory predictions in dark 
energy model with variable Wdc [see Eq. ([37}]. The symbols and lines respectively represent the results with and without the 
EdS approximation: (wojWa) = (—0.9,-0.6) (solid lines and sumbol (—0.9,-1-0.6) (dashed lines and symbol 'x'). In 

bottom panels, the differences between these resuhs are plotted as the fractional error, {P''^'^^'(fc) - p(""™'(fc)}/p(""™^(fc). 
The thick (thin) lines are the results at z — 3(0.5). The vertical arrows indicate the maximum wave number below which the 
prediction of SPT is expected to agree well with the N-body simulations within 1% (see Eq. (I39|l . and Ref. 16]). The filled and 
open arrows are the maximum wave numbers at 2: = 3 and z = 0.5, respectively. 



and [li] 

/ a'' + a« \ 

Wdc{a) = wqWi — — ^ . (38) 

XWiqI + woal ) 

Comparing the numerical results of closure equations with the analytical calculations, we discuss the validity of EdS 
approximation. 

Fig. m shows the one-loop spectra Pii(fc) (le.ft) and P22{k) (right) a.t z — 0.5 and 3, for dark energy model with 
slowly varying Wdo [Eq. (|37p ]. The model parameters wq and Wa were appropriately chosen within the currently 
constrained values of |1 -|- wo\ < 0.1 and \wa\ ^ 0.6 (e.g., Ref. Q). In upper panels, we plot the ratio of power 
spectra, Pab{k) / P^f^(k) , while in lower panels, we plot the fractional difference between the results with and without 
EdS approximation, i.e., {P(i^<^S)(fc) -P(™'^)(fc)}/P(""'°)(fc), where p(EdS) ^nd p(n™) are respectively obtained from 
the analytic and numerical calculations. Similarly, in Fig. O we plot the results in the dark energy model psp . in 
which the equation-of-state parameter has a sharp transition from wi to wq at the scale factor a ~ for a large q. 

The resultant power spectra with EdS approximation underestimate the numerical results without EdS approxima- 
tion in both the density and velocity-divergence part of auto-power spectra. As decreasing the redshift, the deviation 
from numerical results becomes significant, but a level of discrepancy is not so large. These are consistent with the 
previous findings by Refs. (2^. [30^. from the analysis of matter power spectrum. In Figs. [4] and [5l the vertical arrows 
indicate the maximum wave number below which the precision level of one-loop SPT is expected to be better than 
1%. According to Ref. [16), this is determined by solving the following equation with respect to the wavenumber k: 

^ Pti{q;z)dq^ 0.18. (39) 

Note that the maximum wavenumbers given above have been empirically derived by comparison between N-body 
simulations and theoretical predictions, and it seems rather conservative estimates compared to those previously 
proposed [HI, [IS, ISH . Keeping the limitation of the one-loop SPT in mind, we confirm that the analytical treatment 
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FIG. 5: Same as in Fig. 14] but in the dark energy model with equation-of-state parameter psp . The model parameters are 
set to (wo, TOi, fls, (?) = (—1.8, —0.4,0.5,3.41) for set I, and ( — 1.8, —0.4,0.5,25.0) for set II. Note that the case with parameters 
of set //is regarded as extreme one, in which the effective equation-of-state parameter, Wctt = —1 — {2/3)d\n H/dr, sharply 
changes its sign at a ~ 0.5. As for the vertical arrows, see the caption of Fig. ID 



with EdS approximation is a quite good description of the one-loop power spectra and the accuracy of this treatment 
can reach a sub-percent level. This is even true for the model p8|) with the extreme parameter set, i.e., (wq, wi, as, q) — 
(—1.8, —0.4, 0.5, 25.0), in which the effective equation-of-state parameter Wcs = — 1 — {2/3)d\nH/dT, rather than 
Wdc, sharply changes its sign at a ~ 0.5 and eventually approaches WcS ~ 1. 

Therefore, as long as the dark energy models in general relatively are concerned, the analytical calculation with EdS 
approximation is very accurate treatment within the validity range of predictions, and it can give a fast computation 
of the weakly non-linear power spectrum. 



2. Modified gravity models 

Now let us consider the weakly non-linear evolution of the power spectrum in modified gravity models with linear 
Poisson equation. Unlike the dark energy models, the Newton constant is effectively modified, and even the linear 
growth rate generically depends on the scale and time. Thus, the EdS approximation cannot be applied in general 
and the analytical treatment is no longer possible^. Here, we demonstrate that with the use of the present formalism 
and numerical method, the one-loop power spectrum can be accurately computed even in the analytically intractable 
cases. 

We examine two representative modified gravity models whose effective Newton constant manifestly depends on 
the scale and time. One is a phenomenological model in which the Yukawa interaction is added by hand to the 
inverse-square law (e.g., Refs. ^]). The effective Newton constant in this model is given by 



In the DGP model as one of the successful models that explains the late-time cosmic acceleration [^, the effective Newton constant 
Geff depends only on time at the linear-order level, and the analytical calculation of one-loop spectrum is possible with a help of EdS 
approximation, even in the presence of non-linearity in Poisson equation 
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FIG. 6: Weakly non-linear power spectrum at z = 0.5 in a modified gravity model with Yukawa-type interaction (left) and f{R) 
gravity model (right). The symbols, '+', 'x' and '*' respectively represent the resultant spectra in the linear theory, and lines 
the non-linear numerical solutions. The upper panels show the power spectra divided by the linear ones, and the lower ones 
indicate the deviations from the spectra calculated in the flat ACDM model. Note that the symbols in lower panels represent 
the linear theory predictions. 



The parameter A is the characteristic (proper) length at which the Newton force is modified, and the amplitude a 
represents the strength of the deviation from the inverse-square law on large scales. Note that cosmological constraints 
on these parameters have been obtained recently from the galaxy power spectrum of the Sloan Digital Sky Survey 
[33I. [33|. Based on this, we compute the power spectra for specific parameters with (A, a) — (20/i^^Mpc, 1) and 
(100/i-iMpc, 1). 

As another example, we consider the f{R) gravity model. This model has been recently attracted as one of the 
successful models that explains late-time cosmic acceleration [H, [H, [13] (and see also Ref. [1| for a review). The 
f{R) gravity model is given by the generalization of the Einstein-Hilbert action to include arbitrary function of the 
scalar curvature R: 



R + f{R) 
WttG 



(41) 



with Cm being the Lagrangian of ordinary matter. Under the quasi-static treatment relevant for the scales well-inside 
the Hubble horizon, the effective Newton constant becomes (e.g., Ref. [ssj ) 



Geff(fc,r) =G 



1 



3 {k/af+jf 



(42) 



where the quantity jf' is the effective mass of the new scalar degree of freedom, /a, and is defined hy = {(1 + 
fR)/{dfR/dR) — i?}/3. Note that the barred quantity Ji^ impHes the one evaluated in terms of the background 
quantities. In general, the corrections coming from non-linear interaction terms appear in Gcff, but we do not 
consider here. In the present paper, we specifically consider the function /(i?) of the form, f{R) cx R/{AR+ 1) 
[3^, [H, [3^, . In the cosmologically interesting setup with i? and f{R) — > 0, this can be expanded as 



/(i?)~-16^GpA-/flo ( ^ 



(43) 
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The energy density pA is related with the constant A, and Rq and fnQ are the background curvature and the field 
value given by fniRo) at present time. Here, we consider the cases with \fRo\ <^ 1, in which the last term at the 
right-hand side of Eq. ([43|) is safely negligible and the background expansion just follows the same expansion history 
as in the ACDM model. 

Fig.[6]shows the numerical results of one-loop power spectra given dX z — 0.5. Left and right panels plot the results 
for the models with effective Newton constant pO]) and (|42|) , respectively. The upper panels show the ratio of matter 
power spectrum, Pii{k)/ Pl{{k), while in lower panels, the fractional enhancement relative to the ACDM model, i.e., 
{P(fc) — PACDuik)} / PACDuik) , is plotted. In model with Eq. the modification of the gravitational- force law 

appears on large scales, and the effective Newton constant becomes GeS — > (l + a) G. On the other hand, the Newton 
constant on small scales becomes 4/3 times greater than that on large scales in the model with Eq. (|42|l . This scale- 
dependent nature qualitatively explains the results seen in the lower panels, and because of this, the resultant shape 
of the one-loop spectra is significantly altered. Nevertheless, when normalized by the linear power spectra, which 
intrinsically possesses the scale-dependent nature of GeS through the linear growth rate, the differences in the mode 
transfer efficiency between two models turn out to be small (see upper panel). This indicates that the modification 
of gravitational-force law imprinted in the linear power spectrum can be preserved in the weakly non-linear regime, 
and the linear growth rate becomes an important clue to distinguish between various modified gravity models. This 
would be even true for a large class of the modified gravity models with non-linear Poisson equation. 

Finally, it is interesting to note that in the model with Eq. (j42|), there appears the crossing point at which the 
dependence of the ratio P{k)/ P^{k) on |/i^ol is changed. As shown in the upper- right panel, the ratio decreases with 
I //JO I on large scales, while it eventually increases on small scales. This behavior basically reflects the fact that the 
strong gravity on small scales efficiently promotes the mode transfer from the low-fc to high-fc modes. 

VI. DISCUSSION AND CONCLUSION 

In thispaper, on the basis of the non-perturbative framework of the cosmological perturbation theory developed 
by Ref . [l3| , we have presented a numerical scheme to treat non-linear evolution of matter power spectrum. The gov- 
erning equations for matter power spectra are a closed set of evolution equations coupled with non-linear propagator, 
which has been previously derived by truncating the infinite chain of moment equations, with a help of perturbative 
calculation called closure approximation. The present formulation is equivalent to the one-loop level of renormalized 
perturbation theory, and the non-pcrturbative effects of gravitational clustering are effectively incorporated into the 
solution of closure equations. Note that the closure equations consistently reproduces the so-called one- loop results 
of standard perturbation theory if we replace the quantities in the non-linear terms with linear-order ones. The 
numerical scheme presented here can be used for the predictions of matter power spectra in both quasi non-linear and 
non-linear regimes, and is applicable to the analytically intractable cases. The modification to the gravity sector is 
straightforward. 

We have demonstrated that the full non-linear treatment of the closure equations has a ability to treat non-linear 
evolution of power spectrum beyond the validity regime of previous analytical calculations. The resultant shape of the 
non-linear spectrum resembles the N-body result, and the agreement between these two results reaches the accuracy 
of ~ 4%(8%) level at z = 3(z = 0.5). We then focused on the perturbative treatment of closure equations, and 
presented the numerical results of one-loop SPT in various situations. We discussed the validity of the analytical 
treatment based on the Einstein-de Sitter approximation which has been frequently used in the literature. In the 
dark energy models with two representative equation-of-state parameters ([37]) and (|38p . wc found that the analytical 
calculation with Einstein-de Sitter approximation provides an excellent description for the density and velocity- 
divergence components of the one-loop power spectrum. Within the validity range of one-loop spectra, the accuracy 
of this treatment reaches at a sub-percent level. Also, we have studied the one-loop power spectra in a class of modified 
gravity models, in which the effective Newton constant manifestly depends on the scale and time, and the analytical 
calculation is no longer possible. We demonstrated that the scale-dependent modification of the gravitational-force 
law alters the power spectrum significantly, but the efficiency of the mode transfer arising from the non-linear mode 
coupling changes only moderately. In this respect, the modification of the gravity imprinted in the linear power 
spectrum would be preserved in the weakly non-linear regime, and the (scale-dependent) linear growth rate may be 
an important clue to distinguish between various modified gravity models. 

The numerical scheme presented here is a first step toward precisely modeling the non-linear evolution of matter 
power spectrum in various situations. Recently, the non-linear spectrum including the massive neutrinos has been 
investigated by Ref. [4lj based on the approach similar to our formalism Incorporating the effect of massive 

neutrinos into the present formalism is rather straightforward, and the closure equations may be used for a non- 
perturbative calculation of matter power spectrum beyond the free-streaming scales. As another direction, one may 
consider the extension of the present formulation to deal with a wide class of modified gravity models with non-linear 
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Poisson equation. Ref. [l^ presented a general formalism to treat such models and explicitly calculated the one-loop 
power spectrum in DGP and f{R) gravity models from the closure equations. The results for full non-linear treatment 
are left for future work, and will be reported elsewhere. 
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APPENDIX A: DETAILS OF NUMERICAL INTEGRATIONS IN EQS.dlT]) AND ([52]) 

In this appendix, we discuss the technical details on the numerical integrations of the non-linear terms in closure 
equations. In the numerical algorithm presented in Sec. IIVI we must evaluate Eqs. (j31[) and (|32p in advance to the 
time evolution. To compute these integrals, the expressions are first rewritten with the form of the two-dimensional 
integral with a help of the symmetry in the integrands. Then, we introduce the elliptic coordinate used in Ref. 0], 
and perform the integration by the trapezoidal rule taking carefully account of the domain of integration. 

The elliptic coordinate is defined as 



^ ^ / sinh ^ sin /i cos < 

k' = — + q, q = - sinhCsin/isin0 I , (Al) 
y cosh C COS fi 



where the vector k is set to be aligned to the third axis. We introduce 

A' = coshC, F = cos^, (A2) 

where X > 1 and —1<Y<1. In Fig. [71 we schematically plot the elliptic coordinate. In left plot, the origin of the 
elliptic coordinate is O, and the elliptic contour represents a surface of C =const, which is mapped to X = X shown 
in right plot. The vector k points from Fi to F2, which correspond to the two foci of the ellipse. Thus an arbitrary 
vector k' and k — k' can be represented by the vector q. 

In the elliptic coordinate, the vertex functions defined in Eq. ([9]) are recasted as 

7ii2(fc-fc ,fc ) = 7ii2(fc-fc,fc) = , (A3) 

1 — XY X^ + — 2 

7i2i(fc-fc',fcO = _Yy ' 7i2i(fc'-fc,fc) ^ 2(A-y)2 ' '^^^^ 

7222 (fc - fe', k') = '^(^,^V'^P , 7222 (fc' - fc, fc) = (f^)' • (A5) 

From now on, we take (X, F, (/)) as the integration variables instead of (k'^, k'y, k'^). Hence the volume element d^k' 
and the integration domain are changed as 

(■■■) = ^ dX dY{---), (A6) 



Note that, since the integrands of Eqs. PT|) and ([5^ are axially symmetric, we can integrate over the azimuthal angle 
4>, yielding a factor 27r. The lower and upper limits of the integration are determined from the domains of definition 
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of Eqs. and which gives 



Yi{X) =niiiJ -X 



YoiX) = nmx<^ -X 



2fcm — 1 



X 



2k. 



n+l 



Xi — max < 1 



km+l + k; 



■n+l 



(A7) 
(A8) 
(A9) 
(AlO) 



where fc_i and fcAr+i are assigned to ko — fc,nin and fc^v = fcmaxj respectively. With the above preparation, the three- 
dimensional integration pTjl and are reduced to two-dimensional integrations over the domain (|A7|) - (jA10[) . An 
example of integration domain is depicted as a shaded deficient rectangle in right plot of Fig. [71 As mentioned in 
Sec. IIVI we implement the trapezoidal rule to integrate Eqs. ([51]) and (|32p in this domain. 



Xq = max < 1 , 





C = const. 



FIG. 7: The geometrical relation between k' (or q) space and {X, Y) coordinates. The ellipic contour representing a surface 
with ^ = const, in left plot is mapped to X = X in right plot. Circled numbers in both plots corresponds to each other, 
and filled areas, too. Note that the origin of fe'-space is Fi, while that of the elliptic coordinate is O. The striped deficient 
rectangular is an actual domain defined by Eqs. HA7|) - (|A10|I . 



APPENDIX B: CONVERGENCE TEST 



In this appendix, we check the convergence of numerical results obtained with the numerical scheme mentioned in 
Sec. IIVI The test calculations have been done in the ACDM model by varying some numerical parameters. Particularly 
we focus on the initial time of the time evolution, Zinit, and the cutoff wave- number on small scales, fcmax, introduced 
in (|25p . which are most sensitive parameters to the final results. 

The upper panel in Fig. [8] shows the fractional errors between the matter power spectrum Pii(k) for /cmax = 
1,2,5/iMpc"^ and that for fc,nax = lOft-Mpc"^ denoted by P^f , that is, (Pn - P^f)/Plf. In these calculations, 
the logarithmic interval A(logfc) is fixed. In the lower panel, we show the fractional errors for the calculations with 
Zinit = 50 — 300 from the one with Zjnit = 400. Also in these calculations, we fixed the time interval. At. The arrows 
on the horizontal axis are given by Eq. (j39|) . 

Both plots indicate the good convergence of the numerical results in the sense that the fractional errors from each 
reference become smaller as fcmax and ^Tinit increase. We found that we can keep the fractional error sufficiently smaller 
than 1% as long as fcmax is larger than 5/iMpc~^, and Zmit is larger than 200. Particularly, as for ^mit, if we take later 
time, the resultant power spectrum at low-redshifts is harmed by the fact that we neglect the decaying modes in the 
initial conditions [see Eq. (|33p ]. 
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Considering the above results, we fixed kj^ax = 5/i Mpc~ and Zinit = 200 for all numerical calculations presented 
in this paper. Additionally the number of time steps, Nt-, and the wave number bins, Nf^, are chosen as iV^ — 172 
and Nk — 200, respectively, so that the fractional errors are suppressed to a sub-percent level. Moreover, for the 
integration (|A6p . we use a 200 x 200 discrete grid on the integration domain defined by Eqs. (jA7p - (|A10p . 
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FIG. 8: The dependence of the cutoff wave number, fcmax, and the initial time, zinit, on the matter power spectrum. Upper panel 
: The fractional errors of the power spectrum for fcmax = 5^Mpc^^ (solid), femax = 2/iMpc~^ (dashed) and fcmax = l/iMpc~^ 
(dotted) by reference to fcmax ~ lO^Mpc"^. Lower panel : The fractional errors for 2;init = 300 (solid), Zinit = 200 (dashed), 
Zinit = 100 (dotted) and Zinit = 50 (dot-dashed) by reference to zinit ~ 400. As for the vertical arrows, see the caption of Fig.|4] 
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